
*******************************************************************************
*** Program to Produce Table 1 & 2 of Appendix with Bootstrap Standard Errors *
*******************************************************************************

clear
drop _all
matrix drop _all
set more off

cd "$US_Ineq_Repl"
local Processed = "$US_Ineq_Repl/Processed"
local FFL = "$US_Ineq_Repl/Scripts/FFL"

// customize your file path
cd "`FFL'"
global data "`c(pwd)'\Data"
global table "`c(pwd)'\Tables"

local f="$data\"+"9010temptab1.dta"
capture confirm file "`f'"
if c(rc) {
	di "File `f' does not exist. Startiterations from 1"
	local j=1
}
else{
	
	foreach stat in 9010 5010 9050 gini ginilog {
		local f="$data\"+"`stat'temptab1.dta"
		use `f', replace
		count
		local j3`stat' = `r(N)'+1
		mkmat * , matrix(Rt1)
		matrix Rt1`stat'=(nullmat(Rt1`stat')\Rt1)
		
		local f="$data\"+"`stat'temptab2.dta"
		use `f', replace
		count
		local j4`stat' = `r(N)'+1
		mkmat * , matrix(Rt2)
		matrix Rt2`stat'=(nullmat(Rt2`stat')\Rt2)
	}
	
	local j3 = .
	local j4 = .

	foreach stat in 9010 5010 9050 gini ginilog {
		if `j3`stat'' < `j3' | missing(`j3') {
			local j3 = `j3`stat''
		}
		if `j4`stat'' < `j4' | missing(`j4') {
			local j4 = `j4`stat''
		}
	}
	if (`j3'<`j4'){
		local j=`j3'
	}
	else{
		local j=`j4'
	}
	di "current itteration `j'"
}

*** import data ***

global expvar union pubsect manuf nonwhite female partte married smsa ///
	exper exper2 expe2 expe3 expe4 ///
    educ2 educ3 educ4 educ5 edex ee2 ee3 ee4 ee5 ee6 ee7 ee9 ee10 ee11 ee13 ee14 ee15 ///
    indu2 indu3 indu4 indu5 indu6 indu7 indu8 indu9 indu10 indu11 indu12 indu13 indu14 ///
	indu15 indu16 indu17 indu18 indu19 indu20 ///
	state2 state3 state4 state5 state6 state7 state8 state9 ///
	state10 state11 state12 state13 state14 state15 state16 ///
	state17  state18 state19 state20 state21 state22 state23 ///
	state24 state25 state26 state27 state28 state29 state30 ///
	state31 state32 state33 state34 state35 state36 state37 ///
	state38 state39 state40 state42 state43 state44 state45 ///
	state46 state47 state48 state49 state50 state51

use $data\us86.dta , clear
append using $data\us15.dta
append using $data\us86rw15.dta
keep lrwage3 time $expvar orgwgt	// use imputed wages if top-coded
compress
order lrwage3 $expvar
save $data\temp012.dta , replace

use $data\temp012.dta , clear
gen rif_10=.
gen rif_50=.
gen rif_90=.
gen rif_9010=.
gen rif_5010=.
gen rif_9050=.
gen rif_gini=.
gen rif_ginilog=.

*** RIF-decomposition with bootstrap ***
*ssc install rif 
local reps=1000

while `j'<=`reps' {
  preserve
	
  if `j'>1 {
	quietly: bsample, strata(time)  // use original sample in 1st loop
  }

  * calculate RIF for inequality measures
	
  quietly {
  
    // quantile gaps
	forvalues t=0/2 {
	  forvalues qt=10(40)90 {
		local qc = `qt'/100
		pctile valx=lrwage3 [aweight=orgwgt] if time==`t', nq(100)  // 100 percentile grids
		kdensity lrwage3 [aweight=orgwgt] if time==`t', at(valx) ///
		  gen(evalt denst) width(0.065) nograph  // find density at desired grid
		replace rif_`qt'=evalt[`qt']+`qc'/denst[`qt'] if lrwage3>=evalt[`qt'] & time==`t'
		replace rif_`qt'=evalt[`qt']-(1-`qc')/denst[`qt'] if lrwage3<evalt[`qt'] & time==`t'
		drop valx evalt denst
	  }
    }
	replace rif_9010=rif_90-rif_10
	replace rif_5010=rif_50-rif_10
	replace rif_9050=rif_90-rif_50
	
	// variance and gini
	gen wage_gini=exp(lrwage3)
	gen wage_ginilog=lrwage3
	forvalues t=0/2 {
	  foreach stat in gini ginilog {
	  	if ("`stat'"=="ginilog"){
			local detail = "gini"
		}
		else{
			local detail = "`stat'"
		}
		local flag=0
		local c=0
		while `flag'!=1{
			capture rifhdreg wage_`stat' $expvar [aweight=orgwgt] if time==`t', ///
			  rif(`detail') retain(rif`t'_`stat')
			if "`e(cmd)'"=="rifhdreg" {
				local flag=1
			}
			else{
				quietly: bsample, strata(time)
				local c = `c' + 1
				display `c'
				if `c'==3{
					local flag=1
				}
			}
		}
		*rifreg wage_`stat' $expvar [aweight=orgwgt] if time==`t', ///
	    *  `stat' retain(rif`t'_`stat')
	    replace rif_`stat'=rif`t'_`stat' if time==`t'
		drop rif`t'_`stat'
	  }
    }
	replace rif_gini=rif_gini*100
	replace rif_ginilog=rif_ginilog*100
	
  }  // end of quietly

  * OB decomposition
  foreach stat in 9010 5010 9050 gini ginilog {
	display  "doing rep `j', stat `stat'"
	
	// decomposition without reweighing [E(X_1|t=1)- E(X_0|t=0)]B_0
	local flag=0
	while `flag'!=1 {
		quietly: capture oaxaca rif_`stat' $expvar [aweight=orgwgt] if time==0 | time==1, ///
			detail(groupun:union, ///
			groupoth: nonwhite female  married smsa exper-expe4 edex-ee15 state2-state51, /// 
			grouped: educ2-educ5, groupocc: pubsect manuf partte, ///
			groupind: indu2-indu20) weight(0) swap by(time) relax
		
		if "`e(cmd)'"=="oaxaca" {
			local flag=1
		}
		else{
			quietly: bsample, strata(time)
		}
	}
	matrix Ra=e(b)
	di "oaxaca 1/3 completed..."
	
	// composition effects with reweighing [E(X_0|t=1)- E(X_0|t=0)]B_c
	local flag=0
	while `flag'!=1 {
		quietly: capture oaxaca rif_`stat' $expvar [aweight=orgwgt] if time==0 | time==2, ///
			detail(groupun:union, ///
			groupoth: nonwhite female  married smsa exper-expe4 edex-ee15 state2-state51, /// 
			grouped: educ2-educ5, groupocc: pubsect manuf partte, ///
			groupind: indu2-indu20) weight(0) swap by(time) relax
		
		if "`e(cmd)'"=="oaxaca" {
			local flag=1
		}
		else{
			quietly: bsample, strata(time)
		}
	}
	
	matrix Rc=e(b)
	di "oaxaca 2/3 completed..."
	
	// wage structure effects E(X_1|t=1)*[B_1-B_c]
	local flag=0
	while `flag'!=1 {
		quietly: capture oaxaca rif_`stat' $expvar [aweight=orgwgt] if time==1 | time==2, ///
			detail(groupun:union, ///
			groupoth: nonwhite female  married smsa exper-expe4 edex-ee15 state2-state51, /// 
			grouped: educ2-educ5, groupocc: pubsect manuf partte, ///
			groupind: indu2-indu20) weight(0) swap by(time) relax
		
		if "`e(cmd)'"=="oaxaca" {
			local flag=1
		}
		else{
			quietly: bsample, strata(time)
		}
	}

	matrix Rw=-1*e(b)
	di "oaxaca 3/3 completed..."
	
	* collect results

	// w/o reweighting
	matrix Rt1=Ra[1,3..16],Ra[1,6]+Ra[1,11],Ra[1,7]+Ra[1,12],Ra[1,8]+Ra[1,13], ///
		Ra[1,9]+Ra[1,14],Ra[1,10]+Ra[1,15]
	matrix Rt1`stat'=(nullmat(Rt1`stat')\Rt1)

	// w/ reweighting
	matrix Rt2=Ra[1,3],Rc[1,4],Rw[1,5],Rc[1,6..10],Rc[1,5],Rw[1,11..16],Rw[1,4], ///
		Rc[1,6]+Rw[1,11],Rc[1,7]+Rw[1,12],Rc[1,8]+Rw[1,13],Rc[1,9]+Rw[1,14], ///
		Rc[1,10]+Rw[1,15]
	matrix Rt2`stat'=(nullmat(Rt2`stat')\Rt2)
	
	// Saving every 5 iterations
	if (mod(`j',5)==0){
		di "Saving `stat'"
		clear
		svmat double Rt1`stat', name(b) 
		local f="$data\"+"`stat'"+"temptab1.dta"
		save "`f'" , replace
		clear
		svmat double Rt2`stat', name(b) 
		local f="$data\"+"`stat'"+"temptab2.dta"
		save "`f'" , replace
	}
  }  // end of inequality measures loop
  restore
  local j=`j'+1
}  // end of bootstrap loop


*** Get matrices from files and total reps

matrix drop _all
clear

foreach stat in 9010 5010 9050 gini ginilog {
	local f="$data\"+"`stat'temptab1.dta"
	*use if _n <=100 using `f'
	use `f', replace
	count
	local j3`stat' = `r(N)'
	
	mkmat * , matrix(Rt1)
	matrix Rt1`stat'=(nullmat(Rt1`stat')\Rt1)
	matrix drop Rt1
	
	local f="$data\"+"`stat'temptab2.dta"
	*use if _n <=100 using `f'
	use `f', replace
	count
	local j4`stat' = `r(N)'
	
	mkmat * , matrix(Rt2)
	matrix Rt2`stat'=(nullmat(Rt2`stat')\Rt2)
	matrix drop Rt2
}

local j3 = .
local j4 = .

foreach stat in 9010 5010 9050 gini ginilog {
	if `j3`stat'' < `j3' | missing(`j3') {
		local j3 = `j3`stat''
	}
	if `j4`stat'' < `j4' | missing(`j4') {
		local j4 = `j4`stat''
	}
}
if (`j3'<`j4'){
	local reps=`j3'
}
else{
	local reps=`j4'
}

*** define program to export and format results ***

capture program drop mat2table_bs
program define mat2table_bs
  local matnum `1'  // matrix to export: 1 or 2
  local reps `2'   // number of repetitions for bootstrap
	
  *local matnum=3
  *local reps=100
  foreach stat in 9010 5010 9050 gini ginilog {
	clear
	svmat double Rt`matnum'`stat', name(b)  //export coefficients

	* calculate bootstrap standard errors
	local j=1  // loop over rows of output table
	
	foreach var of varlist b* {
		sum `var'
		gen se`j'=`r(sd)' in 1
		replace se`j'=sqrt(`reps'/(`reps'-1))*se`j' in 1  // adj. degree of freedom
		order se`j', after(`var')  // put s.e right "below" coef
		local j=`j'+1
		}
	keep in 1  // use coef from 1st loop (original sample)
	keep b* se*
	
	xpose, clear varname  //transpose row to column
	rename _varname varname
	
	* format results

	// round coef and s.e. by 0.001
	gen y`stat'=string(round(v1,0.001))
	
	// add leading zeros
	replace y`stat'=regexr(y`stat',"^[.]","0.")
	replace y`stat'=regexr(y`stat',"^-[.]","-0.")

	// add parentheses to s.e.
	replace y`stat'="("+y`stat'+")" if substr(varname,1,2)=="se"
	
	// calculate p-values and add aesthetics next to coef
	gen pval=2*normal(-abs(v1/v1[_n+1])) if substr(varname,1,1)=="b"
	replace y`stat'=y`stat'+"*" if pval<=0.1
	replace y`stat'=y`stat'+"*" if pval<=0.05
	replace y`stat'=y`stat'+"*" if pval<=0.01
	keep varname y`stat'
		
	// accumulate results (columns of table)
	merge 1:1 _n using $table\temptablebs.dta , nogen
	sleep 1000 
	save $table\temptablebs.dta , replace	
  }
end

*** format and export table 1 ***
  
clear  // empty file to store results
save $table\temptablebs.dta , emptyok replace
clear
mat2table_bs 1 `reps'
order varname y9010 y5010 y9050 ygini yginilog

// format column titles
label var varname "Variables"
label var y9010 "90-10"
label var y5010 "50-10"
label var y9050 "90-50"
label var ygini "Gini"
label var yginilog "GiniLog"

// format variable labels
replace varname="" if substr(varname,1,2)=="se"
replace varname="Total Change" if varname=="b1"
replace varname="Composition" if varname=="b2"
replace varname="Wage Structure" if varname=="b3"
replace varname="X:Union" if varname=="b4"
replace varname="X:Other" if varname=="b5"
replace varname="X:Education" if varname=="b6"
replace varname="X:Occupation" if varname=="b7"
replace varname="X:Industry" if varname=="b8"
replace varname="W:Union" if varname=="b9"
replace varname="W:Other" if varname=="b10"
replace varname="W:Education" if varname=="b11"
replace varname="W:Occupation" if varname=="b12"
replace varname="W:Industry" if varname=="b13"
replace varname="W:Constant" if varname=="b14"
replace varname="T:Union" if varname=="b15"
replace varname="T:Other" if varname=="b16"
replace varname="T:Education" if varname=="b17"
replace varname="T:Occupation" if varname=="b18"
replace varname="T:Industry" if varname=="b19"

// export table
export excel using $table\tab1bs.xls , replace first(varlabel)
capture rm $table\temptablebs.dta

*** format and export table 2 ***

clear  // empty file to store results
save $table\temptablebs.dta , emptyok replace
clear
mat2table_bs 2 `reps'
order varname y9010 y5010 y9050 ygini yginilog

// format column titles
label var varname "Variables"
label var y9010 "90-10"
label var y5010 "50-10"
label var y9050 "90-50"
label var ygini "Gini"
label var yginilog "GiniLog"

// format variable labels
replace varname="" if substr(varname,1,2)=="se"
replace varname="Total Change" if varname=="b1"
replace varname="Composition" if varname=="b2"
replace varname="Wage Structure" if varname=="b3"
replace varname="X:Union" if varname=="b4"
replace varname="X:Other" if varname=="b5"
replace varname="X:Education" if varname=="b6"
replace varname="X:Occupation" if varname=="b7"
replace varname="X:Industry" if varname=="b8"
replace varname="X:Specification Error" if varname=="b9"
replace varname="W:Union" if varname=="b10"
replace varname="W:Other" if varname=="b11"
replace varname="W:Education" if varname=="b12"
replace varname="W:Occupation" if varname=="b13"
replace varname="W:Industry" if varname=="b14"
replace varname="W:Constant" if varname=="b15"
replace varname="W:Reweighting Error" if varname=="b16"
replace varname="T:Union" if varname=="b17"
replace varname="T:Other" if varname=="b18"
replace varname="T:Education" if varname=="b19"
replace varname="T:Occupation" if varname=="b20"
replace varname="T:Industry" if varname=="b21"

// export table
export excel using $table\tab2bs.xls , replace first(varlabel)

sleep 1000 
capture rm $table\temptablebs.dta


